      implicit real*8 (a-h,o-z)
c
c This program was used for all the simulation experiments in
c James G. MacKinnon and Matthew D. Webb, "Pitfalls when estimating
c treatment effects with clustered data", in The Political Methodologist.
c
c Copyright 2017 James G. MacKinnon
c
      real*8 xmat(10000,51), yvect(10000), ymean(10000)
      real*8 btrue(51), beta(51), betar(51), xpy(51)
      real*8 xpxi(51,51), xpxir(51,51), bpval(4), exparg(200)
      real*8 resid(10000), crvmat(51,51), rres(10000)
      real*8 probs(3,2), crits(3,2), rej(3,7)
      integer year(10000), gsize(200), irej(3,7)
      integer*8 iseed, jseed, isboot
      character*18 mname(7)
      nomax = 10000
      nvmax = 51
      open(unit=8, file='pitfalls.out')
      write(6,*) 'Enter number of replications: '
      read(5,*) nrep
      onrep = 1.d0/nrep
      write(8,100) nrep
 100  format('Number of replications = ', i7)
      write(6,*) 'Enter number of groups and number of observations',
     &  ' per group:'
      read(5,*) ngrp, nperg
      nobs = ngrp*nperg
      write(8,101) ngrp, nperg, nobs
 101  format('Groups: ', i4, ' Obs. per group: ', i5, ' Total number'
     &  ' of observations: ', i6)
      write(6,*) 'Enter gamma parameter for group sizes: '
      read(5,*) gamma
      write(8,*) ' '
      write(8,102) gamma
 102  format('Gamma parameter = ', f8.2)
      if (gamma.eq.0.d0) then
        do ig=1,ngrp
          gsize(ig) = nperg
        end do
        mings = nperg
        maxgs = nperg
      else
        mings = 100000
        maxgs = 0
        rngrp = ngrp
        etot = 0.d0
        do j=1,ngrp
          arg = dfloat(j)/rngrp
          exparg(j) = dexp(gamma*arg)
          etot = etot + exparg(j)
        end do
        mgsum = 0
        ngrp1 = ngrp - 1
        do j=1,ngrp1
          exparg(j) = exparg(j)/etot
          gsize(j) = nobs*exparg(j)
          mgsum = mgsum + gsize(j)
          if (gsize(j).lt.mings) mings = gsize(j)
          if (gsize(j).gt.maxgs) maxgs = gsize(j)
        end do
        gsize(ngrp) = nobs - mgsum
        if (gsize(ngrp).lt.mings) mings = gsize(ngrp)
        if (gsize(ngrp).gt.maxgs) maxgs = gsize(ngrp)
        write(8,103) ngrp, mings, maxgs
 103    format('Sizes of ', i3, ' groups. Minimum: ', i4,
     &    ' Maximum: ', i5)
        write(8,104) (gsize(j),j=1,ngrp)
 104    format(15i5)
      end if
      write(6,*) 'Enter within-group correlation for errors:'
      read(5,*) rhog
      if (rhog.lt.1.d0) then
        sigv = sqrt(rhog)
        sigeps = sqrt(1.d0 - rhog)
      else
        sigv = 1.d0
        sigeps = 0.d0
      end if
      write(8,109) rhog, sigeps, sigv
 109  format('Within-group correlation: ', f7.4,
     &  ' sigeps and sigv: ', 2f8.5) 
      write(6 ,*) 'Enter three 64-bit seeds: '
      read(5,*) iseed, jseed, isboot
      write(8,108) iseed, jseed
 108  format('Main seeds: ', 2i14)
      ntest = 3
      write(6,*) 'How many bootstraps?'
      read(5,*) nboot
      if (nboot.gt.0) then
        ntest = 7
        write(8,112) nboot, isboot
 112    format('Bootstrap samples: ', i4, ' seed: ', i14)
        write(6,*) 'WCR (1), WCU (2), WR (3), WU (4), PC (5)'
        read(5,*) imeth
      end if
      write(6,*) 'Pure treatment model (1) or DiD model (2)?'
      read(5,*) ipure
      if (ipure.eq.1) then
        nvar = 2
        btrue(1) = 1.d0
        btrue(2) = 0.d0
        do i=1,nobs
          xmat(i,1) = 1.d0
          xmat(i,2) = 0.d0
          ymean(i) = btrue(1)
        end do
      else
        write(6,*) 'Enter number of years + first and last to treat:'
        read(5,*) nyear, iyfirst, iylast
        write(8,113) nyear
 113    format('Number of years for DiD model =', i4)
        write(8,114) iyfirst, iylast
 114    format('First and last years to start treatment are', 2i4)
        iynum = iylast + 1 - iyfirst
        nvar = nyear + 1
        i = 0
        do ig=1,ngrp
          jj = 0
          do ii=1,gsize(ig)
            i = i + 1
            jj = jj + 1
            if (jj.gt.nyear) jj = 1
            year(i) = jj
          end do
        end do
        do j=1,nyear
          btrue(j) = 1.d0
        end do
        do i=1,nobs
          ymean(i) = btrue(1)
          do j=1,nyear
            xmat(i,j) = 0.d0
          end do
          xmat(i,year(i)) = 1.d0
        end do
      end if
      nvar1 = nvar - 1
      probs(1,1) = 0.995d0
      probs(2,1) = 0.975d0
      probs(3,1) = 0.95d0
      probs(1,2) = 0.99d0
      probs(2,2) = 0.95d0
      probs(3,2) = 0.90d0
      mname(1) = 'CV1 t-stat Sym'
      mname(2) = 'CV1 t-stat UT'
      mname(3) = 'CV1 t-stat LT'
      write(6,*) 'Enter starting G1, ending G1, increment:'
      read(5,*) igstart, igend, iginc
      do ntrt=igstart,igend,iginc
        if (ipure.eq.1) then
          i = 0
          do ig=1,ntrt
            do ii=1,gsize(ig)
              i = i + 1
              xmat(i,nvar) = 1.d0
            end do
          end do
        end if
        do j=1,3
          do ktest=1,ntest
            irej(j,ktest) = 0
          end do
        end do
        ndfx = ngrp - 1
        rndfx = dfloat(ndfx)
        do j=1,3
          do k=1,2
            call tcdfrinv(probs(j,k),crits(j,k),rndfx)
          end do
        end do
        call timer(it0)
        i1 = 0
        itrep = 1
        do irep=1,nrep
          call timer(it1)
          itx = it1 - it0
          if (itx.gt.1000) then
            time = .01d0*itx
            itrep = irep - itrep
            write(6,66) irep, itrep, time
 66         format(' Replication ', i7, '. Time for last ',i6,
     &        ' replications = ', f8.2, ' seconds.')
            itrep = irep
            it0 = it1
          end if
c
c If DiD model, choose years to treat and create treatment regressor
c
          if (ipure.gt.1) then
            do i=1,nobs
              xmat(i,nvar) = 0.d0
            end do
            i = 0
            do ig=1,ntrt
              call pick64(iseed,iynum,jj)
              iystart = iyfirst + jj - 1
              do ii=1,gsize(ig)
                i = i + 1
                if (year(i).ge.iystart) xmat(i,nvar) = 1.d0
              end do
            end do
          end if
c
c generate yvect
c
          call geny(nobs,nomax,ngrp,iseed,jseed,gsize,sigeps,rhog,
     &      ymean,yvect)
c
c Compute OLS and CRVE
c
          call olsqc(nobs,nomax,nvar,nvmax,i1,ssr,beta,xpy,xpxi,xmat,
     &      yvect,resid)
          itype = 2
          call crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xpxi,resid,
     &      gsize,crvmat)
          tstat = beta(nvar)/sqrt(crvmat(nvar,nvar))
          do j=1,3
            if (tstat.gt.crits(j,1).or.tstat.lt.-crits(j,1)) then
              irej(j,1) = irej(j,1) + 1
            end if
            if (tstat.gt.crits(j,2)) then
              irej(j,2) = irej(j,2) + 1
            end if
            if (tstat.lt.-crits(j,2)) then
              irej(j,3) = irej(j,3) + 1
            end if
          end do
c
c Bootstrap using whatever method was requested
c
          if (nboot.gt.0) then
            if (imeth.eq.1) then
              mname(4) = 'WCR Rad. Sym'
              mname(5) = 'WCR Rad. ET'
              mname(6) = 'WCR Rad. UT'
              mname(7) = 'WCR Rad. LT'
            else if (imeth.eq.2) then
              mname(4) = 'WCU Rad. Sym'
              mname(5) = 'WCU Rad. ET'
              mname(6) = 'WCU Rad. UT'
              mname(7) = 'WCU Rad. LT'
            else if (imeth.eq.3) then
              mname(4) = 'WR Rad. Sym'
              mname(5) = 'WR Rad. ET'
              mname(6) = 'WR Rad. UT'
              mname(7) = 'WR Rad. LT'
            else if (imeth.eq.4) then
              mname(4) = 'WU Rad. Sym'
              mname(5) = 'WU Rad. ET'
              mname(6) = 'WU Rad. UT'
              mname(7) = 'WU Rad. LT'
            else if (imeth.eq.5) then
              mname(4) = 'Pairs Clust. Sym'
              mname(5) = 'Pairs Clust. ET'
              mname(6) = 'Pairs Clust. UT'
              mname(7) = 'Pairs Clust. LT'
            end if
            nwght = 2
            if (imeth.eq.1) then
              bzero = 0.d0
              call olsqc(nobs,nomax,nvar1,nvmax,i1,ssr,betar,xpy,
     &          xpxir,xmat,yvect,rres)
              call wcrboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,
     &          isboot,gsize,yvect,rres,xmat,xpxi,tstat,bzero,bpval)
            else if (imeth.eq.2) then
              bnull = beta(nvar)
              call wcuboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,
     &          isboot,gsize,yvect,resid,xmat,xpxi,tstat,bnull,bpval)
            else if (imeth.eq.3) then
              bzero = 0.d0
              call olsqc(nobs,nomax,nvar1,nvmax,i1,ssr,betar,xpy,
     &          xpxir,xmat,yvect,rres)
              call wrboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,
     &          isboot,gsize,yvect,rres,xmat,xpxi,tstat,bzero,bpval)
            else if (imeth.eq.4) then
              bnull = beta(nvar)
              call wuboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,
     &          isboot,gsize,yvect,resid,xmat,xpxi,tstat,bnull,bpval)
            else if (imeth.eq.5) then
              bnull = beta(nvar)
              call pairsboot(nomax,nvar,nvmax,ngrp,nboot,isboot,
     &          gsize,yvect,xmat,tstat,bnull,bpval)
            end if
            do j=1,3
              prej = 2.d0*(1.d0 - probs(j,1))
              if (bpval(1).lt.prej) irej(j,6) = irej(j,6) + 1
              if (bpval(2).lt.prej) irej(j,7) = irej(j,7) + 1
              if (bpval(3).lt.prej) irej(j,5) = irej(j,5) + 1
              if (bpval(4).lt.prej) irej(j,4) = irej(j,4) + 1
            end do
c
c Cumulate rejections in irej
c
          end if
        end do
        write(8,*) ' '
        do ktest=1,ntest
          do j=1,3
            rej(j,ktest) = irej(j,ktest)*onrep
          end do
          write(8,60) mname(ktest), rhog, ntrt, ngrp,
     &      (rej(j,ktest),j=1,3)
 60       format(a18, f7.4, 2i4, 3f9.6)
        end do
      end do
      stop
      end
      subroutine geny(nobs,nomax,ngrp,iseed,jseed,gsize,
     &  sigeps,rhog,ymean,yvect)
c
c Generate the regressand.
c
c Copyright 2017 James G. MacKinnon
c
      implicit real*8 (a-h,o-z)
      integer*8 iseed, jseed
      real*8 yvect(nomax), uvect(nomax), vvect(ngrp)
      real*8 ymean(nomax), error(nomax)
      integer gsize(ngrp)
      call gg64(iseed,jseed,ngrp,vvect)
      call gg64(iseed,jseed,nobs,uvect)
      i = 0
      sigv = sqrt(rhog)
      do ig=1,ngrp
        do jj=1,gsize(ig)
          i = i + 1
          error(i) = sigeps*uvect(i) + sigv*vvect(ig)
        end do
      end do
      do i=1,nobs
        yvect(i) = ymean(i) + error(i)
      end do
      return
      end      
      subroutine crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xtxinv,
     &  resid,gsize,crvmat)
      implicit real*8 (a-h,o-z)
      real*8 xmat(nomax,nvmax), xtxinv(nvmax,nvmax), resid(nomax)
      real*8 crvmat(nvmax,nvmax), midmat(nvmax,nvmax)
      real*8 xpu(nvmax)
      integer gsize(ngrp)
c
c This routine constructs CRVE covariance matrices for OLS, assuming that
c data are ordered by cluster.
c
c arguments:
c nobs = number of observations
c nvar = number of rhs variables
c ngrp = number of groups
c xtxinv = inverse of x'x matrix (assumed available to calling routine)
c resid = vector of residuals
c gsize = integer vector with number of observations for each group
c itype = 0, 1, or 2 for no blowup or two different blowup factors
c crvmat: output
c
c Copyright 2017 James G. MacKinnon
c
c Preliminary stuff
c
      do j=1,nvar
        do k=j,nvar
          midmat(j,k) = 0.d0
        end do
      end do
      ratio1 = 1.d0
      rngrp = ngrp
      if (itype.eq.1) then
        ratio1 = rngrp/(rngrp - 1.d0)
      else if (itype.eq.2) then
        rnobs = nobs
        ratio1 = (rngrp*(rnobs-1.d0))/((rngrp - 1.d0)*(rnobs-nvar))
      end if
c
c find the middle matrix
c
      igstart = 1
      do ig=1,ngrp
        igend = gsize(ig) + igstart - 1
        do k=1,nvar
          xpu(k) = 0.d0
        end do
        do k=1,nvar
          do i=igstart,igend
            xpu(k) = xpu(k) + xmat(i,k)*resid(i)
          end do
        end do
        do j=1,nvar
          do k=j,nvar
            midmat(j,k) = midmat(j,k) + xpu(j)*xpu(k)
          end do
        end do
        igstart = igend + 1
      end do
      do j=1,nvar
        do k=j,nvar
          midmat(k,j) = midmat(j,k)
        end do
      end do
c
c put everything together
c
      do j=1,nvar
        do k=j,nvar
          crvmat(j,k) = 0.d0
          do  jk=1,nvar
            do jl=1,nvar
              add = xtxinv(j,jk)*midmat(jk,jl)*xtxinv(jl,k)
              crvmat(j,k) = crvmat(j,k) + add
            end do
          end do
          crvmat(j,k) = ratio1*crvmat(j,k)
          crvmat(k,j) = crvmat(j,k)
        end do
      end do
      return
      end
      subroutine wcrboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,isboot,
     &  gsize,yvect,rres,xmat,xpxi,tstat,best,bpval)
      implicit real*8 (a-h,o-z)
c
c Calculates WCR wild cluster bootstrap P value for hypothesis that
c the last coefficient is zero.
c Use various weights, based on nwght.
c nwght = 1 is normal
c nwght = 2 is Rademacher
c nwght = 3 is uniform
c nwght = 4 is Mammen
c nwght = 6 is Webb
c
c Copyright 2017 James G. MacKinnon
c
      real*8 yvect(nomax), xmat(nomax,nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
      real*8 crvmat(nvmax,nvmax), vboot(nomax), rres(nomax)
      real*8 bboot(nvmax), bpval(4)
      real*8 wild(6), wp(6), yrfit(nomax), yboot(nomax), resboot(nomax)
      integer gsize(ngrp)
      integer*8 isboot, jsboot
      itype = 2
      if (nwght.eq.1) then
        jsboot = isboot/2
      end if
      if (nwght.eq.3) then
        root12 = sqrt(12.d0)
      end if
      if (nwght.eq.6) then
        wild(1) = sqrt(1.5d0)
        wild(2) = 1.d0
        wild(3) = sqrt(0.5d0)
        wild(4) = - wild(3)
        wild(5) = -1.d0
        wild(6) = - wild(1)
        wprob = 1.d0/6.d0
        wp(1) = wprob
        wp(2) = 2.d0*wprob
        wp(3) = 3.d0*wprob
        wp(4) = 4.d0*wprob
        wp(5) = 5.d0*wprob
      end if
      if (nwght.eq.2) then
        wild(1) = 1.d0
        wild(2) = -1.d0
        wprob = 0.5d0
      end if
      if (nwght.eq.4) then
        wprob = (sqrt(5.d0) + 1.d0)/(2.d0*sqrt(5.d0))
        wild(1) = -(sqrt(5.d0) - 1.d0)/2.d0
        wild(2) = 1.d0 - wild(1)
      end if
c
c Use rres to compute yrfit
c
      do i=1,nobs
        yrfit(i) = yvect(i) - rres(i)
      end do
c
c Start the bootstrap loop
c
      bpval(1) = 0.d0
      bpval(2) = 0.d0
      bpval(4) = 0.d0
      atstat = abs(tstat)
      do iboot=1,nboot
c
c Generate bootstrap samples
c
        if (nwght.ne.1) then
          call ggun64(isboot,ngrp,0.d0,1.d0,vboot)
        else
          call gg64(isboot,jsboot,ngrp,vboot)
        end if
        igstart = 1
        do ig=1,ngrp
          igend = gsize(ig) + igstart - 1
c
c Webb six-point weights
c
          if (nwght.eq.6) then
            if (vboot(ig).lt.wp(1)) erand = wild(1)
            if (vboot(ig).ge.wp(5)) erand = wild(6)
            do k=2,5
              if (vboot(ig).lt.wp(k).and.vboot(ig).ge.wp(k-1)) then
                erand = wild(k)
              end if
            end do
          end if
c
c Rademacher weights
c
          if (nwght.eq.2) then
            if (vboot(ig).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
          end if
c
c Mammen weights
c
          if (nwght.eq.4) then
            if (vboot(ig).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
          end if
c
c Standard normal weights
c
          if (nwght.eq.1) then
            erand = vboot(ig)
          end if
c
c Uniform weights
c
          if (nwght.eq.3) then
            erand = root12*(vboot(ig) - 0.5d0)
          end if
c
          do i=igstart,igend
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
          igstart = igend + 1
        end do
c
c Estimate unrestricted model using bootstrapped data.
c
        i1 = 1
        call olsqc(nobs,nomax,nvar,nvmax,i1,ssr,bboot,xpy,xpxi,xmat,
     &    yboot,resboot)
        call crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xpxi,
     &    resboot,gsize,crvmat)
        btstat = (bboot(nvar) - best)/sqrt(crvmat(nvar,nvar))
        if (btstat.gt.tstat) then
          bpval(1) = bpval(1) + 1.d0
        else
          bpval(2) = bpval(2) + 1.d0
        end if
        abtstat = abs(btstat)
        if (abtstat.gt.atstat) then
          bpval(4) = bpval(4) + 1.d0
        end if
      end do
      bpval(1) = bpval(1)/nboot
      bpval(2) = bpval(2)/nboot
      bpval(3)= 2.d0*min(bpval(1),bpval(2))
      bpval(4) = bpval(4)/nboot
      return
      end
      subroutine wcuboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,isboot,
     &  gsize,yvect,ures,xmat,xpxi,tstat,best,bpval)
      implicit real*8 (a-h,o-z)
c
c Calculates WCU wild cluster bootstrap P value for hypothesis that
c the last coefficient is zero.
c Use various weights, based on nwght.
c nwght = 1 is normal
c nwght = 2 is Rademacher
c nwght = 3 is uniform
c nwght = 4 is Mammen
c nwght = 6 is Webb
c
c Copyright 2017 James G. MacKinnon
c
      real*8 yvect(nomax), xmat(nomax,nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
      real*8 crvmat(nvmax,nvmax), vboot(nomax), ures(nomax)
      real*8 bboot(nvmax), bpval(4)
      real*8 wild(6), wp(6), yfit(nomax), yboot(nomax), resboot(nomax)
      integer gsize(ngrp)
      integer*8 isboot, jsboot
      itype = 2
      if (nwght.eq.1) then
        jsboot = isboot/2
      end if
      if (nwght.eq.3) then
        root12 = sqrt(12.d0)
      end if
      if (nwght.eq.6) then
        wild(1) = sqrt(1.5d0)
        wild(2) = 1.d0
        wild(3) = sqrt(0.5d0)
        wild(4) = - wild(3)
        wild(5) = -1.d0
        wild(6) = - wild(1)
        wprob = 1.d0/6.d0
        wp(1) = wprob
        wp(2) = 2.d0*wprob
        wp(3) = 3.d0*wprob
        wp(4) = 4.d0*wprob
        wp(5) = 5.d0*wprob
      end if
      if (nwght.eq.2) then
        wild(1) = 1.d0
        wild(2) = -1.d0
        wprob = 0.5d0
      end if
      if (nwght.eq.4) then
        wprob = (sqrt(5.d0) + 1.d0)/(2.d0*sqrt(5.d0))
        wild(1) = -(sqrt(5.d0) - 1.d0)/2.d0
        wild(2) = 1.d0 - wild(1)
      end if
c
c Use ures to compute yfit
c
      do i=1,nobs
        yfit(i) = yvect(i) - ures(i)
      end do
c
c Start the bootstrap loop
c
      bpval(1) = 0.d0
      bpval(2) = 0.d0
      bpval(4) = 0.d0
      atstat = abs(tstat)
      do iboot=1,nboot
c
c Generate bootstrap samples
c
        if (nwght.ne.1) then
          call ggun64(isboot,ngrp,0.d0,1.d0,vboot)
        else
          call gg64(isboot,jsboot,ngrp,vboot)
        end if
        igstart = 1
        do ig=1,ngrp
          igend = gsize(ig) + igstart - 1
c
c Webb six-point weights
c
          if (nwght.eq.6) then
            if (vboot(ig).lt.wp(1)) erand = wild(1)
            if (vboot(ig).ge.wp(5)) erand = wild(6)
            do k=2,5
              if (vboot(ig).lt.wp(k).and.vboot(ig).ge.wp(k-1)) then
                erand = wild(k)
              end if
            end do
          end if
c
c Rademacher weights
c
          if (nwght.eq.2) then
            if (vboot(ig).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
          end if
c
c Mammen weights
c
          if (nwght.eq.4) then
            if (vboot(ig).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
          end if
c
c Standard normal weights
c
          if (nwght.eq.1) then
            erand = vboot(ig)
          end if
c
c Uniform weights
c
          if (nwght.eq.3) then
            erand = root12*(vboot(ig) - 0.5d0)
          end if
c
          do i=igstart,igend
            yboot(i) = yfit(i) + ures(i)*erand
          end do
          igstart = igend + 1
        end do
c
c Estimate unrestricted model using bootstrapped data.
c
        i1 = 1
        call olsqc(nobs,nomax,nvar,nvmax,i1,ssr,bboot,xpy,xpxi,xmat,
     &    yboot,resboot)
        call crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xpxi,
     &    resboot,gsize,crvmat)
        btstat = (bboot(nvar) - best)/sqrt(crvmat(nvar,nvar))
        if (btstat.gt.tstat) then
          bpval(1) = bpval(1) + 1.d0
        else
          bpval(2) = bpval(2) + 1.d0
        end if
        abtstat = abs(btstat)
        if (abtstat.gt.atstat) then
          bpval(4) = bpval(4) + 1.d0
        end if
      end do
      bpval(1) = bpval(1)/nboot
      bpval(2) = bpval(2)/nboot
      bpval(3)= 2.d0*min(bpval(1),bpval(2))
      bpval(4) = bpval(4)/nboot
      return
      end
      subroutine wrboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,isboot,
     &  gsize,yvect,rres,xmat,xpxi,tstat,best,bpval)
      implicit real*8 (a-h,o-z)
c
c Calculates WR wild bootstrap P value for hypothesis that the last
c coefficient is zero using CV1.
c Use various weights, based on nwght.
c nwght = 1 is normal
c nwght = 2 is Rademacher
c nwght = 3 is uniform
c nwght = 4 is Mammen
c nwght = 6 is Webb
c
c Copyright 2017 James G. MacKinnon
c
      real*8 yvect(nomax), xmat(nomax,nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
      real*8 crvmat(nvmax,nvmax), vboot(nomax), rres(nomax)
      real*8 bboot(nvmax), bpval(4)
      real*8 wild(6), wp(6), yrfit(nomax), yboot(nomax), resboot(nomax)
      integer gsize(ngrp)
      integer*8 isboot, jsboot
      itype = 2
      if (nwght.eq.1) then
        jsboot = isboot/2
      end if
      if (nwght.eq.3) then
        root12 = sqrt(12.d0)
      end if
      if (nwght.eq.6) then
        wild(1) = sqrt(1.5d0)
        wild(2) = 1.d0
        wild(3) = sqrt(0.5d0)
        wild(4) = - wild(3)
        wild(5) = -1.d0
        wild(6) = - wild(1)
        wprob = 1.d0/6.d0
        wp(1) = wprob
        wp(2) = 2.d0*wprob
        wp(3) = 3.d0*wprob
        wp(4) = 4.d0*wprob
        wp(5) = 5.d0*wprob
      end if
      if (nwght.eq.2) then
        wild(1) = 1.d0
        wild(2) = -1.d0
        wprob = 0.5d0
      end if
      if (nwght.eq.4) then
        wprob = (sqrt(5.d0) + 1.d0)/(2.d0*sqrt(5.d0))
        wild(1) = -(sqrt(5.d0) - 1.d0)/2.d0
        wild(2) = 1.d0 - wild(1)
      end if
c
c Use rres to compute yrfit
c
      do i=1,nobs
        yrfit(i) = yvect(i) - rres(i)
      end do
c
c Start the bootstrap loop
c
      bpval(1) = 0.d0
      bpval(2) = 0.d0
      bpval(4) = 0.d0
      atstat = abs(tstat)
      do iboot=1,nboot
c
c Generate bootstrap samples
c
        if (nwght.ne.1) then
          call ggun64(isboot,nobs,0.d0,1.d0,vboot)
        else
          call gg64(isboot,jsboot,nobs,vboot)
        end if
c
c Webb six-point weights
c
        if (nwght.eq.6) then
          do i=1,nobs
            if (vboot(i).lt.wp(1)) erand = wild(1)
            if (vboot(i).ge.wp(5)) erand = wild(6)
            do k=2,5
              if (vboot(i).lt.wp(k).and.vboot(i).ge.wp(k-1)) then
                erand = wild(k)
              end if
            end do
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
c
c Rademacher weights
c
        else if (nwght.eq.2) then
          do i=1,nobs
            if (vboot(i).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
c
c Mammen weights
c
        else if (nwght.eq.4) then
          do i=1,nobs
            if (vboot(i).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
c
c Standard normal weights
c
        else if (nwght.eq.1) then
          do i=1,nobs
            erand = vboot(i)
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
c
c Uniform weights
c
        else if (nwght.eq.3) then
          do i=1,nobs
            erand = root12*(vboot(i) - 0.5d0)
            yboot(i) = yrfit(i) + rres(i)*erand
          end do
        end if
c
c Estimate unrestricted model using bootstrapped data.
c
        i1 = 1
        call olsqc(nobs,nomax,nvar,nvmax,i1,ssr,bboot,xpy,xpxi,xmat,
     &    yboot,resboot)
        call crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xpxi,
     &    resboot,gsize,crvmat)
        btstat = (bboot(nvar) - best)/sqrt(crvmat(nvar,nvar))
        if (btstat.gt.tstat) then
          bpval(1) = bpval(1) + 1.d0
        else
          bpval(2) = bpval(2) + 1.d0
        end if
        abtstat = abs(btstat)
        if (abtstat.gt.atstat) then
          bpval(4) = bpval(4) + 1.d0
        end if
      end do
      bpval(1) = bpval(1)/nboot
      bpval(2) = bpval(2)/nboot
      bpval(3)= 2.d0*min(bpval(1),bpval(2))
      bpval(4) = bpval(4)/nboot
      return
      end
      subroutine wuboot(nobs,nomax,nvar,nvmax,ngrp,nwght,nboot,isboot,
     &  gsize,yvect,ures,xmat,xpxi,tstat,best,bpval)
      implicit real*8 (a-h,o-z)
c
c Calculates WU wild bootstrap P value for hypothesis that the last
c coefficient is zero using CV1.
c Use various weights, based on nwght.
c nwght = 1 is normal
c nwght = 2 is Rademacher
c nwght = 3 is uniform
c nwght = 4 is Mammen
c nwght = 6 is Webb
c
c Copyright 2017 James G. MacKinnon
c
      real*8 yvect(nomax), xmat(nomax,nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
      real*8 crvmat(nvmax,nvmax), vboot(nomax), ures(nomax)
      real*8 bboot(nvmax), bpval(4)
      real*8 wild(6), wp(6), yfit(nomax), yboot(nomax), resboot(nomax)
      integer gsize(ngrp)
      integer*8 isboot, jsboot
      itype = 2
      if (nwght.eq.1) then
        jsboot = isboot/2
      end if
      if (nwght.eq.3) then
        root12 = sqrt(12.d0)
      end if
      if (nwght.eq.6) then
        wild(1) = sqrt(1.5d0)
        wild(2) = 1.d0
        wild(3) = sqrt(0.5d0)
        wild(4) = - wild(3)
        wild(5) = -1.d0
        wild(6) = - wild(1)
        wprob = 1.d0/6.d0
        wp(1) = wprob
        wp(2) = 2.d0*wprob
        wp(3) = 3.d0*wprob
        wp(4) = 4.d0*wprob
        wp(5) = 5.d0*wprob
      end if
      if (nwght.eq.2) then
        wild(1) = 1.d0
        wild(2) = -1.d0
        wprob = 0.5d0
      end if
      if (nwght.eq.4) then
        wprob = (sqrt(5.d0) + 1.d0)/(2.d0*sqrt(5.d0))
        wild(1) = -(sqrt(5.d0) - 1.d0)/2.d0
        wild(2) = 1.d0 - wild(1)
      end if
c
c Use ures to compute yfit
c
      do i=1,nobs
        yfit(i) = yvect(i) - ures(i)
      end do
c
c Start the bootstrap loop
c
      bpval(1) = 0.d0
      bpval(2) = 0.d0
      bpval(4) = 0.d0
      atstat = abs(tstat)
      do iboot=1,nboot
c
c Generate bootstrap samples
c
        if (nwght.ne.1) then
          call ggun64(isboot,nobs,0.d0,1.d0,vboot)
        else
          call gg64(isboot,jsboot,nobs,vboot)
        end if
c
c Webb six-point weights
c
        if (nwght.eq.6) then
          do i=1,nobs
            if (vboot(i).lt.wp(1)) erand = wild(1)
            if (vboot(i).ge.wp(5)) erand = wild(6)
            do k=2,5
              if (vboot(i).lt.wp(k).and.vboot(i).ge.wp(k-1)) then
                erand = wild(k)
              end if
            end do
            yboot(i) = yfit(i) + ures(i)*erand
          end do
c
c Rademacher weights
c
        else if (nwght.eq.2) then
          do i=1,nobs
            if (vboot(i).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
            yboot(i) = yfit(i) + ures(i)*erand
          end do
c
c Mammen weights
c
        else if (nwght.eq.4) then
          do i=1,nobs
            if (vboot(i).le.wprob) then
              erand = wild(1)
            else
              erand = wild(2)
            end if
            yboot(i) = yfit(i) + ures(i)*erand
          end do
c
c Standard normal weights
c
        else if (nwght.eq.1) then
          do i=1,nobs
            erand = vboot(i)
            yboot(i) = yfit(i) + ures(i)*erand
          end do
c
c Uniform weights
c
        else if (nwght.eq.3) then
          do i=1,nobs
            erand = root12*(vboot(i) - 0.5d0)
            yboot(i) = yfit(i) + ures(i)*erand
          end do
        end if
c
c Estimate unrestricted model using bootstrapped data.
c
        i1 = 1
        call olsqc(nobs,nomax,nvar,nvmax,i1,ssr,bboot,xpy,xpxi,xmat,
     &    yboot,resboot)
        call crve(nobs,nomax,nvar,nvmax,ngrp,itype,xmat,xpxi,
     &    resboot,gsize,crvmat)
        btstat = (bboot(nvar) - best)/sqrt(crvmat(nvar,nvar))
        if (btstat.gt.tstat) then
          bpval(1) = bpval(1) + 1.d0
        else
          bpval(2) = bpval(2) + 1.d0
        end if
        abtstat = abs(btstat)
        if (abtstat.gt.atstat) then
          bpval(4) = bpval(4) + 1.d0
        end if
      end do
      bpval(1) = bpval(1)/nboot
      bpval(2) = bpval(2)/nboot
      bpval(3)= 2.d0*min(bpval(1),bpval(2))
      bpval(4) = bpval(4)/nboot
      return
      end
      subroutine pairsboot(nomax,nvar,nvmax,ngrp,nboot,isboot,
     &  gsize,yvect,xmat,tstat,best,bpval)
      implicit real*8 (a-h,o-z)
c
c This routine performs the pairs cluster bootstrap for hypothesis that
c the last coefficient is zero using CV1.
c
c It is a bit tricky, because sample sizes vary across bootstrap
c samples, and because the X matrix may not have full rank
c
c Copyright 2017 James G. MacKinnon
c
      real*8 yvect(nomax), xmat(nomax,nvmax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
      real*8 crvmat(nvmax,nvmax)
      real*8 bboot(nvmax), bpval(4)
      real*8, dimension(:,:), allocatable :: xboot
      real*8, dimension(:), allocatable :: yboot, resboot
      integer gsize(ngrp), gsum(ngrp), agsize(ngrp), iosave(ngrp)
      integer*8 isboot
      character*20 estcoef
      itype = 2
      gsum(1) = 0
      do ig=2,ngrp
        gsum(ig) = gsum(ig-1) + gsize(ig-1)
      end do
      bpval(1) = 0.d0
      bpval(2) = 0.d0
      bpval(4) = 0.d0
      atstat = abs(tstat)
      do iboot=1,nboot
c
c Pick the clusters to include in the bootstrap dataset
c
 887    continue
        nobsb = 0
        do igrp=1,ngrp
          call pick64(isboot,ngrp,iout)
          agsize(igrp) = gsize(iout)
          iosave(igrp) = iout
          nobsb = nobsb + agsize(igrp)
        end do
        allocate (yboot(nobsb))
        allocate (resboot(nobsb))
        allocate (xboot(nobsb,nvar))
        nobsb = 0
        do igrp=1,ngrp
          nobsbp = nobsb
          nobsb = nobsb + agsize(igrp)
          istart = gsum(iosave(igrp))
          do ii=1,agsize(igrp)
            i = nobsbp + ii
            yboot(i) = yvect(istart + ii)
            do k=1,nvar
              xboot(i,k) = xmat(istart + ii,k)
            end do
          end do
        end do
        i1 = 0
        call olsqc(nobsb,nobsb,nvar,nvmax,i1,ssr,bboot,xpy,xpxi,xboot,
     &    yboot,resboot)
c
c Check whether the coefficient can be computed
c
        write(estcoef,888) bboot(nvar)
 888    format(d20.14)
        if (estcoef(2:4).eq.'NaN') then
          deallocate (yboot,resboot,xboot)
          go to 887
        end if
        call crve(nobsb,nobsb,nvar,nvmax,ngrp,itype,xboot,xpxi,
     &    resboot,agsize,crvmat)
        std = sqrt(crvmat(nvar,nvar))
        if (std.lt.1.d-10) then
          deallocate (yboot,resboot,xboot)
          go to 887
        end if
        btstat = (bboot(nvar) - best)/std
        if (btstat.gt.tstat) then
          bpval(1) = bpval(1) + 1.d0
        else
          bpval(2) = bpval(2) + 1.d0
        end if
        abtstat = abs(btstat)
        if (abtstat.gt.atstat) then
          bpval(4) = bpval(4) + 1.d0
        end if
        deallocate (yboot,resboot,xboot)
      end do
      bpval(1) = bpval(1)/nboot
      bpval(2) = bpval(2)/nboot
      bpval(3)= 2.d0*min(bpval(1),bpval(2))
      bpval(4) = bpval(4)/nboot
      return
      end
      subroutine olsqc(nobs,nomax,nvar,nvmax,i1,ssr,beta,
     & xpy,xpxi,xmat,yvect,resid)
      implicit real*8 (a-h,o-z)
c
c This Cholesky-based olsq routine is fast, but may not be as
c accurate as it could be.
c
c nomax is dimension of things that logically have length nobs (so that
c   nomax => nobs); nvmax is dimension of things that logically have
c length nvar (so that nvmax => nvar).
c
c This version returns beta, ssr, resid.
c It also returns xpxi, and can be told to use an xpxi given to it.
c It does not compute the covariance matrix; multiply xpxi by the
c square of sigest to do that
c
c Copyright 2002 James G. MacKinnon
c
      real*8 xmat(nomax,nvmax), yvect(nomax), beta(nvmax), resid(nomax)
      real*8 xpy(nvmax), xpxi(nvmax,nvmax)
c
c  form and invert xpx matrix unless already done
c
      if(i1.ne.0) go to 250
c only one triangle of xpxi computed.
      do i=1,nvar
      do j=i,nvar
        sum = 0.d0
        do k=1,nobs
          sum = sum + xmat(k,i)*xmat(k,j)
        end do
        xpxi(j,i) = sum
        xpxi(i,j) = sum
      end do
      end do
c
c  now invert xpx
c
      call chol2(xpxi,nvmax,nvar,kf2)
  250 continue
c
c preliminary work has now been done.
c
c form xpy
c
      do i=1,nvar
        sum = 0.d0
        do k=1,nobs
          sum = sum + xmat(k,i)*yvect(k)
        end do
        xpy(i) = sum
      end do
c
c  now form estimates of beta.
c
      do i=1,nvar
        beta(i) = 0.d0
        do j=1,nvar
          beta(i) = beta(i) + xpxi(i,j)*xpy(j)
        end do
      end do
c
      ssr = 0.d0
      do k=1,nobs
        sum = yvect(k)
        do i=1,nvar
          sum = sum - xmat(k,i)*beta(i)
        end do
        ssr = ssr + sum**2
        resid(k) = sum
      end do
      return
      end
      subroutine chol2(a,m,n,kf2)
c
c This routine uses the cholesky decomposition to invert a real
c symmetric matrix.
c Copyright 2002 James G. MacKinnon
c
      implicit real*8 (a-h,o-z)
      real*8 a(m,m)
      kf2 = 0
      do 8 i=1,n
        kl = i - 1
        do 7 j=i,n
          if (i.gt.1) then
            do k=1,kl
              a(i,j) = a(i,j) - a(k,i)*a(k,j)
            end do
          else
            if (a(i,i).le.0.d0) then
              kf2 = i
              go to 20
            end if
          end if
          if (i.eq.j) then
            a(i,i) = dsqrt(a(i,i))
          else
            if (j.eq.i+1) ooa = 1.d0/a(i,i)
            a(i,j) = a(i,j)*ooa
          end if
 7      continue
 8    continue
      do 13 i=1,n
        do j=i,n
          ooa = 1.d0/a(j,j)
          if (i.ge.j) then
            t = 1.d0
            go to 12
          end if
          kl = j - 1
          t = 0.d0
          do k=i,kl
            t = t - a(i,k)*a(k,j)
          end do
 12       a(i,j) = t*ooa
        end do
 13   continue
      do 16 i=1,n
        do 15 j=i,n
          t = 0.d0
          do k=j,n
            t = t + a(i,k)*a(j,k)
          end do
          a(i,j) = t
          a(j,i) = t
 15     continue
 16   continue
 20   return
      end
      subroutine tcdfrinv(prob,crit,rndf)
c
c This routine inverts the Student's t CDF using bisection
c It allows degrees of freedom to be real can calls tcdfreal.f
c Copyright 2014 James G. MacKinnon
c
      implicit real*8 (a-h,o-z)
      data toler/1.d-10/
      if (prob.gt.0.5d0) then
        aprob = 1.d0 - prob
      else
        aprob = prob
      end if
      top = 0.d0
      ftop = .5d0
      bot = -20.d0
      if (rndf.lt.2.d0) bot = -50.d0
      if (rndf.lt.1.5d0) bot = -100.d0
      if (rndf.lt.1.0d0) bot = -200.d0
      if (rndf.lt.0.80d0) bot = -500.d0
      if (rndf.lt.0.70d0) bot = -10000.d0
      if (rndf.lt.0.50d0) bot = -50000.d0
      if (rndf.lt.0.40d0) bot = -500000.d0
      call tcdfreal(bot,fbot,rndf)
      do i=1,1000
        xmid = .5d0*(bot + top)
        call tcdfreal(xmid,value,rndf)
        fmid = value - aprob
        if (fmid.lt.0.d0) then
          bot = xmid
          fbot = fmid
        else
          top = xmid
          ftop = fmid
        end if
        diff = abs(bot - top)
        if (diff.lt.toler) go to 100
      end do
 100  continue
      crit = .5d0*(bot + top)
      if (aprob.ne.prob) crit = - crit
      return
      end
      subroutine tcdfreal(arg,value,rndf)
c
c This routine finds c.d.f. of Student's t distribution with real D.F.
c arg = argument (double precision)
c value = value of c.d.f. (double precision)
c rndf = degrees of freedom (double precision)
c
c Calls functions betai, betacf, gammln, which are also in tcdf.f,
c so just include tcdf.f.
c
      implicit real*8 (a-h,o-z)
      a = .5d0*rndf
      b = .5d0
      x = rndf/(rndf + arg**2)
      value = .5d0*betai(a,b,x)
      if (arg.ge.0.d0) value = 1.d0 - value
      return
      end
      subroutine tcdf(arg,value,ndf)
c
c This routine finds c.d.f. of Student's t distribution.
c arg = argument (double precision)
c value = value of c.d.f. (double precision)
c ndf = degrees of freedom (integer)
c
c Calls functions betai, betacf, gammln
c
      implicit real*8 (a-h,o-z)
      a = .5d0*ndf
      b = .5d0
      x = ndf/(ndf + arg**2)
      value = .5d0*betai(a,b,x)
      if (arg.ge.0.d0) value = 1.d0 - value
      return
      end
      function betai(a,b,x)
      implicit real*8 (a-h,o-z)
      if (x.lt.0.d0.or.x.gt.1.d0) then
        write(6,*) 'Bad argument x in betai.'
      end if
      if (x.eq.0.d0.or.x.eq.1.d0) then
        bt = 0.d0
      else
        bt = exp(gammln(a+b) - gammln(a) - gammln(b)
     1  + a*dlog(x) + b*dlog(1.d0-x))
      endif
      if (x.lt.(a+1.d0)/(a+b+2.d0)) then
        betai = bt*betacf(a,b,x)/a
        return
      else
        betai = 1.d0 - bt*betacf(b,a,1.d0-x)/b
        return
      endif
      end
      function gammln(xx)
c
c This routine calculates the log of the Gamma function.
c
      implicit real*8 (a-h,o-z)
      real*8 cof(6),stp,half,one,fpf,x,tmp,ser
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     *    -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*log(tmp) - tmp
      ser=one
      do 11 j=1,6
        x = x + one
        ser = ser + cof(j)/x
11    continue
      gammln = tmp + log(stp*ser)
      return
      end
      function betacf(a,b,x)
      implicit real*8 (a-h,o-z)
      parameter (itmax=200,eps=1.d-7)
      am = 1.d0
      bm = 1.d0
      az = 1.d0
      qab = a+b
      qap = a+1.d0
      qam = a-1.d0
      bz = 1.d0 - qab*x/qap
      do 11 m = 1,itmax
        em = m
        tem = em+em
        d = em*(b-m)*x/((qam+tem)*(a+tem))
        ap = az+d*am
        bp = bz+d*bm
        d = -(a+em)*(qab+em)*x/((a+tem)*(qap+tem))
        app = ap+d*az
        bpp = bp+d*bz
        aold = az
        am = ap/bpp
        bm = bp/bpp
        az = app/bpp
        bz = 1.d0
        if (abs(az-aold).lt.eps*abs(az)) then
          betacf = az
          return
        end if
11    continue
      write(6,*) 'In betacf, a or b too big, or itmax too small.'
        betacf = az
      return
      end
      subroutine gg64(iseed,jseed,nobs,uvect)
c
c This routine generates a vector of N(0,1) random variables
c using two 64-bit random number generators using multiplicative
c congruential schemes; transformed to normal by Marsaglia-Bray
c Copyright 2008 James G. MacKinnon
c
      implicit real*8 (a-h,o-z)
      implicit integer*8 (i-n)
      integer i, ii, nr, nobs, nobs2
      real*8 uvect(nobs)
      data ia/2307085864/, im/9223372036854775783/
      data iq/3997845152/, ir/214644455/, ofm/1.084202172485504D-19/
      data ia2/2946716567/, im2/9223372036854775097/
      data iq2/3130050626/, ir2/1671854155/, ofm2/1.084202172485504D-19/
c
c 1 / 9223372036854775783 = 1.084202172485504E-19
c 1 / 9223372036854775097 = 1.084202172485504E-19
c
      nobs2 = nobs/2
      nr = nobs - 2*nobs2
      i = 0
      do 10 ii=1,nobs2
        i = i + 2
c start of generating r1
    1   ihi = iseed/iq
        ilo = iseed - ihi*iq
        iseed = ia*ilo - ir*ihi
        if (iseed.le.0) iseed = iseed + im
        r1 = iseed*ofm
c end of generating r1; start of generating r2
        ihi = jseed/iq2
        ilo = jseed - ihi*iq2
        jseed = ia2*ilo - ir2*ihi
        if (jseed.le.0) jseed = jseed + im2
        r2 = jseed*ofm2
c end of generating r2
        v1 = 2.d0*r1 - 1.d0
        v2 = 2.d0*r2 - 1.d0
        r = v1**2 + v2**2
        if (r.ge.1.d0) go to 1
        fac = sqrt(-2.d0*log(r)/r)
        uvect(i) = v1*fac
        uvect(i-1) = v2*fac
   10 continue
      if (nr.eq.0) then
        return
      else
c this is for when nobs is odd
        i = i + 1
c start of generating r1
    2   ihi = iseed/iq
        ilo = iseed - ihi*iq
        iseed = ia*ilo - ir*ihi
        if (iseed.le.0)  iseed = iseed + im
        r1 = iseed*ofm
c end of generating r1; start of generating r2
        ihi = jseed/iq2
        ilo = jseed - ihi*iq2
        jseed = ia2*ilo - ir2*ihi
        if (jseed.le.0) jseed = jseed + im2
        r2 = jseed*ofm2
c end of generating r2
        v1 = 2.d0*r1 - 1.d0
        v2 = 2.d0*r2 - 1.d0
        r = v1**2 + v2**2
        if (r.ge.1.d0) go to 2
        fac = sqrt(-2.d0*log(r)/r)
        uvect(i) = v2*fac
      end if
      return
      end
      subroutine ggun64(iseed,nobs,amin,amax,uvect)
c
c This routine generates a vector of uniforms (amin,amax) using
c a 64-bit generator, the first one in gg64.
c Copyright 2013 James G. MacKinnon
c
      real*8 amin, amax, adiff
      real*8 uvect(nobs)
      real*8 ofm, rand
      integer*8 ia, im, iq, ir, iseed, ihi, ilo
      integer i, nobs
      data ia/2307085864/, im/9223372036854775783/
      data iq/3997845152/, ir/214644455/, ofm/1.084202172485504D-19/
c
c 1 / 9223372036854775783 = 1.084202172485504E-19
c
      adiff = amax - amin
      do i=1,nobs
        ihi = iseed/iq
        ilo = iseed - ihi*iq
        iseed = ia*ilo - ir*ihi
        if (iseed.le.0) iseed = iseed + im
        rand = iseed*ofm
        uvect(i) = adiff*rand + amin
      end do
      return
      end
      subroutine pick64(iseed,nobs,iout)
c
c This routine returns a random integer between 1 and nobs
c Copyright 2013 James G. MacKinnon
c
      real*8 ofm, unif
      integer*8 ia, im, iq, ir, iseed, ihi, ilo
      integer nobs, iout
      data ia/2307085864/, im/9223372036854775783/
      data iq/3997845152/, ir/214644455/, ofm/1.084202172485504D-19/
c
c 1 / 9223372036854775783 = 1.084202172485504E-19
c
      ihi = iseed/iq
      ilo = iseed - ihi*iq
      iseed = ia*ilo - ir*ihi
      if (iseed.le.0) iseed = iseed + im
      unif = iseed*ofm
      unif = nobs*unif + 1.d0
      iout = unif
      return
      end

